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Abstract 

A general purpose orthotropic elasto-plastic computational constitutive material model has been 
developed to improve predictions of the response of composites subjected to high velocity impact. The 
three-dimensional orthotropic elasto-plastic composite material model is being implemented initially for 
solid elements in LS-DYNA as MAT213. In order to accurately represent the response of a composite, 
experimental stress-strain curves are utilized as input, allowing for a more general material model that can 
be used on a variety of composite applications. The theoretical details are discussed in a companion 
paper. This paper documents the implementation, verification and qualitative validation of the material 
model using the T800-F3900 fiber/resin composite material. 

Introduction 

Composite materials are now beginning to provide uses hitherto reserved for metals in structural 
systems such as airframes and engine containment systems, wraps for repair and rehabilitation and 
ballistic/blast mitigation systems. While material models exist that can be used to simulate the response of 
materials in these demanding structural systems, they are often designed for use with specific 
homogeneous materials such as metals (Refs. 1 and 2), polymers (Ref. 3) and wood (Ref. 4). While 
material models exist to simulate the impact response of composites, most are tailored specifically for a 
particular application and have limitations - purely elastic, no rate sensitivity, implementation for solid 
elements only, and limited damage and failure characterization (Refs. 5 to 10). 

A few examples of different approaches to composite damage modeling include Littell et al. 

(Ref. 11), Matzenmiller et al. (Ref. 12), and others (Refs. 13 to 18). Yen’s comprehensive model (Ref. 19) 
is implemented in LS-DYNA as *MAT_162. Availability of constitutive models for composite materials 
in LS-DYNA is discussed in some length in our companion paper (Ref. 20) and that explains the 
motivation for this study. 
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In the next section we briefly discuss the theoretical development of the constitutive model. This is 
followed by details of the numerical algorithm, verification tests for a composite used in the aerospace 
industry, the use of the composite in a high speed contact example and finally, concluding remarks 
including details of our ongoing and future work. 


Orthotropic 3D Elasto-Plastic Composite Material Model 

The material law that the model is built upon describes the elastic and permanent deformation of the 
composite with full three-dimensional implementation suitable for solid elements. Only the relevant 
equations used in the numerical algorithm are presented here. The theoretical basis and discussions are 
available in a companion paper (Ref. 20). Current development of the model includes the complete rate 
independent elasto-plastic deformation model, with damage and failure to be added later. 

The commonly used Tsai-Wu composite failure model is generalized as a yield surface for the 
plasticity model and is defined as 
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The stress components of the yield function coefficients correspond to the current yield stresses 
associated with a set of experimental stress-strain curves, which vary with the effective plastic strain, thus 
allowing the model to describe different hardening properties in each direction. The non-associative flow 
surface is defined as 



where the Hfj terms are the constant flow rule coefficients. The rest of the details can be found in the 
companion paper (Ref. 20). 

The first six flow rule coefficients are computed directly from the assumed flow rule coefficient 
value, Hii and the plastic Poisson’s ratios. 
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The last three flow rule coefficients (H 44 , H 55 , Hee) are calculated iteratively with the user-defined 
range to best fit the input experimental shear stress versus shear strain curves (Ref. 20). For example, to 
determine H 44 , the following optimization problem is posed using the input curve for the “1-2” shear case. 
Find H 44 to minimize 


/(//44) = X (^22), -(^12). 
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where n is the number of data points on the master curve, {622 ) i* effective stress value from the 
master curve, and (^ 12 ). effective stress value for the shear curve using an assumed value of H 44 . 


Numerical Algorithm 

In this initial implementation, we do not account for rate and temperature dependence, damage 
accumulation nor failure. The following experimental data are needed as LS-DYNA input. 

(1) A total of 12 true stress versus true strain curves under quasi-static test conditions from (a) uniaxial 
tension in 1, 2 and 3-directions, (b) uniaxial compression in 1, 2 and 3-directions, (c) shear in 1-2, 
2-3 and 3-1 planes and (d) off-axis (e.g., 45°) uniaxial tension or compression in 1-2, 2-3 and 3-1 
planes. The tabulated x-y data for each curve is read in and used appropriately. 

(2) The modulus of elasticity, Poisson’s ratio and average plastic Poisson’s ratio obtained from the 
tension and compression tests. 

The effective plastic strain can be written in terms of the experimental stress versus total strain data. 

In the 1 -direction for example, this is written as 




aliisP 


( 4 ) 


where is the experimental compressive true stress in the 1 -direction, Sn is the total true strain in the 
1 -direction, is the elastic modulus in the 1 -direction, sf’j is the true plastic strain in the 1 -direction. 
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8 ^ is the effective plastic strain and h is the current value of the plastic potential function, which is 
equivalent to the effective stress, shown in Equation ( 2 ). The material model uses a typical elastic stress 
update as 


<j"+'=o"+CAf:(E-i'’) 


( 5 ) 


where C is the orthotropic stiffness matrix. At is the time step, £ is the total strain rate and £^ is the 
plastic strain rate. The stiffness matrix is written in terms of the compliance matrix as 
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where E^j are the elastic moduli in the principal material directions, are the elastic shear moduli and 
Vy are the elastic Poisson’s ratios. The yield stresses used to determine the flow rule coefficients are 
summarized into a single vector, corresponding to each of the experimental test curves as 
q^=[0[, 0^2 'J33 ^22 012 023 CT31 ®«-i2 ^^45-23 045-31] • Lastly, the consistency 

condition is written in terms of the rate of yield function as follows 


G 


'23 


/ = — 0 + — q = 0 

where the terms are as defined earlier. Equation ( 7 ) can be expanded as follows by introducing the 
elasto-plastic constitutive law 


( 7 ) 


V da) 5 q dX 


( 8 ) 


where a is written in terms of the constitutive matrix and strain rates and X is the effective plastic strain. 
Solving for the effective plastic strain rate produces the following equation. 
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To numerically solve for the increment in effective plastic strain, a radial return algorithm is used. 
Initially, a perfectly elastic response is assumed. Thus, an elastic predictor is used to compute the stresses 
as 
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With the elastic trial stresses computed, a trial yield function value can be calculated and used to 
determine if the load step is elastic or plastic. It should be noted that the elastic trial stress is used as the 
stress in the next step if the trial yield function value is less than or equal to zero. Otherwise, the stress 
state is beyond the yield surface and must be brought back using a radial return and plastic corrector, AA. . 
If the trial yield function happens to be greater than zero, then A A must be greater than zero. The value of 

AA is determined using a secant iteration, with AA^ = 0 for the first iteration and the second value for the 
increment of the plastic multiplier determined from consistency. 
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where the derivatives of q are taken as zero, meaning that the response is perfectly plastic. The initial 
calculation of the second increment corresponds to perfect plasticity in order to ensure that the stress state 
returns to the interior of the yield surface, resulting in a negative value of the yield function and bounding 
the solution. From each estimate of AA , the corresponding corrected stresses can be computed as well as 
the yield stresses obtained from the input stress versus effective plastic strain (A) data. With these terms 
calculated for both estimates of AA , the yield function value for each can be computed with the third 
estimate evaluated as 


AA3=AA1-/1^^^^ (12) 

f-f' 

The corrected stresses and yield stresses associated with AA^ are then calculated and used to obtain a 
new estimate of the yield function value. Convergence of the secant iteration is determined by the 
following conditions, where /b is the value of the yield function given the effective plastic strain AA^, 
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with the secant iteration continuing if the new yield function value is not less than some defined tolerance. 
Once convergence is met, the plastic multiplier increment to return the stress state to the yield surface is 
known and the stresses can be updated as 


i"*^=a^-C-AX— 

dn 


(14) 


Finally, the yield stresses are updated as well, using the new value of the overall effective plastic 
strain, X, in each input curve to determine the corresponding yield stress level as 


g’‘*^=q[x''+AX^) 


(15) 


This results in anisotropic strain hardening, as a yield stress increase in each direction is initiated with 
plasticity in any direction, but at different levels. A detailed algorithm is presented below. 

Step 1: Preprocessing 

(a) Convert stress/strain input curves to stress versus effective plastic strain using Equation (4). 

(b) Store initial yield stresses in q . 

(c) Calculate and store the elastic moduli from the initial yield stress and strain values. 

(d) Store the three elastic Poisson’s ratio values and flow rule coefficients from input. 

(e) Compute optimal values of the flow rule coefficients so as to match the input curves as closely as 
possible. 

Step 2: Initialization 

Parameters available at start of step: . 


Step 3: Elastic predictor 

Compute and set the yield function coefficients using Equation (la). 

Construct the constitutive matrix using Equation (6). 

Compute elastic trial stresses, , using Equation (10). 

Compute the trial yield function, , using the elastic trial stresses in Equation (1). 

If < 0 , the elastic solution is correct, set AX^ = 0 and go to stress update. Else go to plastic 
corrector. 

Step 4: Plastic corrector 
Set AA,‘ = 0. 

Calculate A A. from Equation (11). 

Loop through secant iteration for specified max number of iterations: 

Compute the new estimate of the stress for each plastic multiplier increment j using 

Equation (14). 
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Calculate the effective plastic strains ) at the next time step using the estimates of the plastic 


multiplier increment; then obtain the strain hardening parameters with the effective plastic strain 
from the curve data (true stress as a function of effective plastic strain). 

3 3 

Compute new plastic multiplier increment, AX . Calculate new estimate for the yield function, / . 


If 


< tol , set AA. = , exit the loop and go to stress update. Else update secant iteration 


parameters using Equation (14) and go to next step of secant iteration. 


Step 5: Stress Update 

Calculate using Equation (14). 


Step 6: Update history variables for plastic work and work hardening parameters (q and A. ) 

Set Xy^J^l = Xy^+ AXj^ . 

Determine new yield stresses, q„+i, using X^^^i and Equation (15). 


Model Verification 

The composite material model was tested and verified using experimental data obtained from 
T800S/3900-2B [P2352W-19] BMS8-276 Rev-H-Unitape fiber/resin unidirectional composite (Ref. 21). 
Toray describes T800S as an intermediate modulus, high tensile strength graphite fiber. The epoxy resin 
system is labeled E3900 where a toughened epoxy is combined with small elastomeric particles to form a 
compliant interface or interleaf between fiber plies to resist impact damage and delamination (Ref. 22). 
Magnified views of the composite are shown in Eigures 1 and 2. 



Figure 1 . — Side view obtained 
using optical microscopy. 



Figure 2. — Logitudinal view obtained 
using a SEM. 
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The source of the input data for the model (both experimental and virtual) is listed in Table 1. In the 
Data column, Experimental refers to experimental data generated at Wichita State University (Ref. 21). 
For cases where actual experimental data were not available, numerical simulations (virtual experiments) 
were conducted to obtain the required stress-strain curves. To conduct the numerical experiments, 
micromechanics analyses were conducted where the properties of the fiber and matrix were used to 
simulate the response of the overall composite. Two methods were used to conduct the numerical 
experiments. First, the NASA GRC developed code MAC-GMC was used (Ref. 23). MAC-GMC utilized 
the Generalized Method of Cells to conduct a fully analytical simulation of the composite response. A 
second approach used a Virtual Testing Software System (Ref. 24) in which finite element models were 
constructed and analyzed where the fiber and matrix were explicitly modeled. For this particular 
unidirectional composite, as indicated in the table, due to the transverse isotropy of the material the 
3-direction response was set equal to the 2-direction response. Likewise, the 1-3 shear response was set 
equal to the 1-2 shear response, and the off-axis response in the 1-3 plane was set equal to the off-axis 
response in the 1-2 plane. More details on the procedures used in the numerical experiments will be given 
in a future report. 

The fiber (transversely isotropic, linear elastic) and the matrix (isotropic, elastic perfectly plastic) 
properties used in the numerical experiments are listed in Table 2. The properties for the latter are not 
publicly available. Hence they were correlated to obtain the values shown in the table. 

The values of the coefficients in the flow rule were determined from the experimental (actual or 
virtual) data and the procedures outlined in Equations (3a), (3b), and (3c). The correlation curves that 
were used to optimize the shear coefficients in the flow law are shown in Figures 3 and 4. In the plots, the 
[90°] transverse tensile curve was chosen as the baseline master curve for the optimization. 


TABLE 1.— GENERATION OE INPUT DATA EOR 
T800-E3900 COMPOSITE 


Curve 

Data 

Tension Test (I -Direction) 

Experimental 

Tension Test (2-Direction) 

MAC-GMC and VTSS 

Tension Test (3-Direction) 

Transverse isotropy 

Compression Test (I -Direction) 

Experimental 

Compression Test (2-Direction) 

MAC-GMC and VTSS 

Compression Test (3-Direction) 

Transverse isotropy 

Pure Shear Test (1-2 Plane) 

Experimental 

Pure Shear Test (2-3 Plane) 

MAC-GMC and VTSS 

Pure Shear Test (1-3 Plane) 

Transverse isotropy 

Off-Axis Test (45°, 1-2 Plane) 

MAC-GMC and VTSS 

Off-Axis Test (45°, 2-3 Plane) 

MAC-GMC and VTSS 

Off-Axis Test (45°, 1-3 Plane) 

Transverse isotropy 


TABEE 2.— T800-E3900 EIBER AND MATRIX PROPERTIES 
(VOLUME ERACTION = 0.54) 


Property 

Eiber 

Matrix 

El, psi 

4(109 

5(109 

E 2 , psi 

2.25(109 


E 3 , psi 

2.25(109 


V12 

0.2 

0.35 

V23 

0.25 


V13 

0.25 


G, psi 

1.5(109 

1.85(109 

ay, psi 


2(109 
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Effective Plastic Stress (psi) Effective Plastic Stress (psi) 


Optimized H44 Value-Generated Curve 



Effective Plastic Strain 

Figure 3. — Comparison of experimental curve with optimized Haa (and Hee) value curve. 


Optimized Hss Value-Generated Curve 



Effective Plastic Strain 

Figure 4. — Comparison of experimental curve with optimized H 55 value curve. 
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After the determination of the required constants, the verification process was executed. In the 
verification test results shown in the figures that follow, the term Experimental (Input) refers to 
experimental data (when available) or to data created via virtual testing, and the term Simulated 
(x-element) refers to the data created using results from finite element simulations where x is the number 
of finite elements in the FE model. 

Schematics for the tension and compression tests are shown in Figures 5 and 6, respectively. The 
simulations for the tension and compression tests in all three directions were executed using different 
mesh sizes ranging from 1 to 64 elements. The 64-element meshes for the tension and compression test 
cases are shown in Figures 7 and 8, respectively. 



E C 

^ 8” 




-^1 

2” B 

I ^2, 3 


t=0.25” 


A t=0.25” 

(a) (b) 

Figure 5. — Schematics for tension test simulation cases (a) 1 -direction (b) 2 and 3-directions. 


2(y) 



■1 (X) 


(a) 


1 (y) 



- 2/3 (X) 


(b) 


Figure 6.— Schematics for compression test simulation cases (a) 1 -direction (b) 2 and 3-directions. 



Figure 7. — The 64-element mesh for tension test cases. 



Figure 8. — The 64-element mesh 
for compression test cases. 
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For both tension and compression tests, nodes on face ABC were constrained in the x-direction and 
the center node was also constrained in the y-direction. A 0.5 in./s displacement in the x-direction was 
applied to all the nodes on the right face, E-D. For the compression tests, care was taken to ensure that the 
model yielded but did not buckle. The simulated and experimental stress-strain curves for the tension and 
compression tests are shown in Figures 9 and 10. In all cases the correlation between the input curves and 
the model predictions was excellent. Note that the transverse tension and compression curves (both the 
numerical experiments and model simulations) were artificially extended well into the nonlinear range. 
The extension was required in order to both fully characterize the model and to explore the capability of 
the model to correctly simulate the nonlinear response of the composite. 



0.00 0.01 0.02 0.03 0.04 0.05 

Strain 

(a) 


0-Degree Compression Test 



(b) 

Figure 9. — Simulated and experimental stress-strain curves for 1 -direction (a) tension and (b) compression. 
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Stress (psi) 


90-Degree Tension Test (2/3-direction) 
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(a) 


90-Degree Compression Test (2/3-direction) 
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Figure 10. — Simulated and experimental stress-strain curves for 2/3-directions (a) tension and (b) compression. 
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Stress (psi) 


A schematic for the pure shear test is shown in Figure 11. Nodes on the bottom surface are restrained 
in the x and y-directions while nodes on the top surface are restrained in the y-direction. A displacement 
controlled loading of 0.5 in./s in the x-direction is applied to all the nodes on the top surface. The 
64-element mesh is shown in Figure 12. The simulated and experimental stress-strain curves in the 
1-2/3-1 plane are shown in Figure 13 and in the 2-3 plane in Figure 14. Once again, the correlation 
between the input curves and the model predictions was excellent. 

The (tension) off-axis test model in the 1-2 plane is shown in Figure 15. Nodes on face ABC were 
constrained in the x-direction and the center node was also constrained in the y-direction. The local 
element coordinate system was rotated 45° in order to simulate the off-axis response. The simulated and 
experimental stress-strain curves are compared in Figure 16, showing very little differences in the results, 
which at least partially demonstrates the ability of the model simulate the response of a composite when it 
is subjected to a multi-axial stress state (in the material axis system). Finally, Figure 17 shows the results 
from a loading-unloading-reloading test case. The model captured the results expected for a pure 
plasticity formulation, including a linear elastic unload with the original modulus, and a resumption of 
plastic flow upon reloading to the original unload point. Capturing the reduction in the unloading 
modulus expected for materials of this type will require the addition of a complementary damage model 
to the formulation. 


2 

(y) 













.25” 



Figure 1 1 .—Schematic for pure shear test in 
1 -2/3-1 plane. 


1 (X) 



i-X 


Figure 12. — The 64-element mesh for pure shear 
test cases. 


1-2/3-1 Shear Test 



Figure 13. — Simulated and experimental stress-strain curves for pure shear in the 1 -2/3-1 plane. 
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Figure 14. — Simulated and experimental stress-strain curves for pure shear in the 2-3 plane. 
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Figure 15.— Schematic for 45° off axis test in 1 -2/3-1 plane. 
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Figure 16.— Simulated and experimental stress-strain curves for 45° off axis test in 1 -2/3-1 plane. 


Off-Axis Tension (1-2 plane, 45 Degree) 


















Input/E 

Simulai 

ixperimental 
ted (1 element) 
ted (4-element) 
ted (16-element) 
ted (64-element) - 




Simulai 

Simulai 




— -Simulai 


NAS A/TM— 20 15-218446 


14 



0.00 0.01 0.02 0.03 0.04 0.05 

Strain 

Figure 17. — Simulated and experimental stress-strain curves for unloading/reloading in the 2-direction. 

In summary, simulations of the curves represented by the 12 input cases as an initial verification 
study shows very good correlation with the input data. We have used boundary conditions that resemble 
the experimental test setup as closely as possible so as to capture the actual test conditions. The use of 
four different meshes with increasing mesh density, shows that the results are mesh independent. 

Structural Test 

To gage the ability of the model to produce reasonable results under more realistic structural loading 
conditions, a structural test case was designed to qualitatively validate the developed model. A simply- 
supported 12 in. long by 2.4 in. wide by 0.6 in. thick plate made of T800-F3900 composite material is 
subjected to an impact from a steel ball. Figure 18 shows the schematic diagram. The nodes on the bottom 
of the left edge are constrained in the x, y and z-directions while nodes on the bottom of the right edge 
were constrained only in the y-direction. In the absence of experimental data, the primary objectives were 
to ensure that the material model (in the absence of damage and failure) would yield realistic results. That 
is, the ball would impact the plate and bounce back, and the peak stresses in the plate would occur around 
the point of contact immediately after the impact. The LS-DYNA FE model is shown in Figure 19. 
References to LS-DYNA keywords can be found in Hallquist (Ref. 25). 

The model has 6,478 nodes and 5,707 elements of which 1,147 8-noded hexahedral elements are used 
to model the plate (master surface) and the remaining 4,560 elements are used in the ball (slave surface). 
The ball is modeled using *MAT_3 and the plate as *MAT_213. *CONTACT_AUTOMATIC_ 

SURF ACE_TO_SURF ACE is used to describe the contact between the ball and plate. Both the static (fs) 
and the dynamic (fd) coefficients of friction are set to 0.1. A 2 percent viscous damping coefficient is 
used. Segment-based contact is used by setting the soft parameter to 2. The steel ball is given a -700 in./s 
initial velocity in the y-direction. The analysis is carried out for 0.1 sec. The results are as expected. The 
exaggerated deformed shape is shown in Eigure 20. The von Mises stress throughout the test duration is 
plotted for the elements closest to the point of contact in Eigure 21. 
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Figure 18— Schematic diagram of the test. 



Figure 19. — LS-DYNA FE model at the start of the Figure 20. Transient shape after impact, 

analysis. 



Figure 21 . — von Mises stress versus time for elements near point of impact. 
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Conclusions 


A generalized three-dimensional composite material model is developed and verified. The preliminary 
results are encouraging as the implemented constitutive model is able to reproduce the experimental stress- 
strain curves. In addition, a simple contact problem involving the composite material yields reasonable 
results. Future work includes the following: adding damage accumulation and failure in the material model, 
supporting temperature and strain rate dependent stress-strain curves, using alternate methods to compute 
the plastic multiplier when the initial yield and flow surface convexity conditions are not met, using a fiber- 
referenced convected system to track the principal material directions and supporting plate/shell elements 
and to validate the developed model using data from several structural tests. 
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